## Figure 4: Effects of Neighborhood Race and Income on Expected Wait Times

## install.packages(c("tidyverse", "fixest"))
## library(tidyverse)
## library(fixest)

## SET WORKING DIRECTORY HERE
## setwd()

## Loading data
## load("dta.RData")

## Figure 4a
## Across-service
summary(race_expected_across <- feols(log_expected_time ~ factor(white_third) |
                                        city^open_month^open_year, data = dta %>% filter(city != "New York"), cluster = "geo"))
race_expected_across_coefs = as.data.frame(race_expected_across$coeftable) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Within-service
summary(race_expected_within <- feols(log_expected_time ~ factor(white_third) |
                                        city_service^open_month^open_year, data = dta %>% filter(city != "New York"),
                                      cluster = "geo"))
race_expected_within_coefs = as.data.frame(race_expected_within$coeftable) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

race_coefs = race_expected_within_coefs %>%
  bind_rows(., race_expected_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "% White Tercile",
       y = "% Change in Expected Wait \n Time vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_rect(colour = "black"),
        legend.position = "bottom") 

## Figure 4b
## Across-service
summary(inc_expected_across <- feols(log_expected_time ~ factor(inc_third) |
                                       city^open_month^open_year, data = dta %>% filter(city != "New York"),
                                     cluster = "geo"))

inc_expected_across_coefs = as.data.frame(inc_expected_across$coeftable) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Within-service
summary(inc_expected_within <- feols(log_expected_time ~ factor(inc_third) |
                                       city_service^open_month^open_year, data = dta %>% filter(city != "New York"),
                                     cluster = "geo"))
inc_expected_within_coefs = as.data.frame(inc_expected_within$coeftable) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

inc_coefs = inc_expected_within_coefs %>%
  bind_rows(., inc_expected_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "Per Capita Income Tercile",
       y = "% Change in Expected Wait \n Time vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_rect(colour = "black"),
        legend.position = "bottom") 
